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The ^-statistic, derived by Jaranowski, Krolak & Schutz (1998), is the optimal (frequentist) 
statistic for the detection of nearly periodic gravitational waves from known neutron stars, in the 
presence of stationary, Gaussian detector noise. The jF-statistic was originally derived for the 
case of a single detector, whose noise spectral density was assumed constant in time, and for a 
single known neutron star. Here we show how the ^-statistic can be straightforwardly generalized 
to the cases of 1) a network of detectors with time-varying noise curves, and 2) a collection of 
known sources (e.g., all known millisecond pulsars within some fixed distance). Fortunately, all 
the important ingredients that go into our generalized ^"-statistics are already calculated in the 
single-source/single-detector searches that are currently implemented, e.g., in the LIGO Software 
Library, so implementation of optimal multi-detector, multi-source searches should require negligible 
additional cost in computational power or software development. This paper also includes an analysis 
of the likely efficacy of a collection-type search, and derives criteria for deciding which candidate 
sources should be included in a collection, if one is trying to maximize the detectability of the 
whole. In particular we show that for sources distributed uniformly in a thin disk, the strongest 
source in the collection should have signal-to-noise-squared ~ 5 times larger than weakest source, for 
an optimized collection. We show that gravitational waves from collection of the few brightest (in 
gravitational waves) neutron stars could perhaps be detected before the single brightest source, but 
that this is far from guaranteed. Once gravitational waves from the few brightest neutron stars have 
been discovered, grouping more distant (individually undetectable) pulsars into collections, and then 
searching for those collections, should be an effective way of measuring the average gravitational- 
wave strengths of those more distant pulsars. 

PACS numbers: 95.55.Ym, 04.80.Nn, 95.75.Pq, 97.60.Gb 



I. INTRODUCTION 

The ^-statistic, as first derived by Jaranowski, Kro- 
lak & Schutz pj (hereinafter referred to as JKS), is the 
optimal frequentist statistic for the detection of nearly 
periodic gravitational waves (GWs) from a known neu- 
tron star. In the original JKS version, the ^"-statistic was 
derived only for the case of a single GW detector (which 
was assumed to have stationary noise characteristics) and 
a single known neutron star (assumed to be emitting 
GWs at the neutron star's rotation frequency and/or at 
twice its rotation freuqency). Here we show how the T- 
statistic can be generalized in a straightforward manner 
to the cases of 1) a network of detectors with time- varying 
noise curves, and 2) an entire collection of known sources. 
Fortunately, all the important ingredients that go into 
the generalized .F-statistic are already calculated in the 
singlc-detector/single-source searches that are currently 
implemented, e.g., in the LIGO Software Library Q, so 
implementation of optimal multi-detector and/or multi- 
source searches should require negligible additional cost 
in software development and computation. 

We note that the problem of optimally combining data 
from different detectors has already been solved for sev- 
eral types of GW searches. For the case of inspiralling 
binaries, we refer the reader to Bose, Pai & Dhurand- 
har Q and to Finn £|; for the case of GW bursts, to 
Sylvestre and for the case of LISA observations of 
galactic, stellar-mass binaries, to both Krolak et al. 
and Rogan & Bose Q. (LISA can effectively be treated 



as a network of three independent GW detectors.) Our 
analysis in Sec. Ill is especially similar to that of Krolak 
et al. [|| and Rogan & Bose Q , since formally the sources 
considered there are equivalent to GW pulsars. Like GW 
pulsars, the stellar-mass binaries visible to LISA are effec- 
tively monochromatic sources that can be characterized 
by four amplitude parameters, in addition to the GW fre- 
quency and the two angles specifying the source position 
on the sky. 

The basic idea of somehow combining the signals 
from many individually-undetectable sources or events, 
in hopes of finding a statistical excess, is also hardly a 
new one. In GW astronomy, a good example is the sug- 
gestion of looking for GW bursts associated with gamma- 
ray bursts by cross-correlating the outputs of LIGO's 
LI and HI detectors over short time windows coincident 
with hundreds of observed gamma-ray bursts @ . But our 
application of this idea to the population of known mil- 
lisecond pulsars appears to be new. We investigate when 
this strategy is likely to be effective and derive useful cri- 
teria for deciding how many and which sources should 
be included in the collection, in order to maximize the 
detectability of that group. 

The plan of this paper is as follows. In Sec. [H] we 
briefly establish notation; we generally try to align our 
notation with that of JKS, to ease comparison with their 
work. In Sec. IIIII we derive the JF-statistic for a network 
of N detectors and a single source. This multi-detector 
^-statistic follows a \ 2 distribution with 4 degrees of 
freedom, exactly as with the single-detector version. We 
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consider the general case where the detectors have corre- 
lated noises, but of course our expressions simplify in the 
case where noises from different detectors are uncorre- 
cted. As a bonus, our results immediately show how to 
appropriately time- weight the data in the (realistic) case 
that the detector noise spectra are slowly time-varying. 
(The appropriate time- weighting for a single detector has 
already been derived by Itoh et al. and is implemented 
in the LIGO Software Library. We give an independent 
derivation, since here it follows trivially.) 

In Sec. II VI we extend the JF-statistic to the case of a col- 
lection of known sources. If there are M known sources 
(each emitting at a single, known frequency), then the 
correct JF-statistic for the entire collection follows a % 2 
distribution with 4M degrees of freedom. (This is true for 
both the single-detector case and for an iV-detector net- 
work.) The most interesting target population is clearly 
(some subset of) the known millisecond pulsars. We con- 
sider two particularly interesting cases: 1) a collection 
of the few very brightest GW pulsars, and 2) a larger 
collection of more distant GW pulsars. We investigate 
the expected gains from both these types of multi-source 
searches, under the reasonable assumption that there ex- 
ists some population of GW pulsars that is uniformly 
spread thoughout the Galactic disk. As a further illus- 
tration of multi-source searches, we estimate the sensi- 
tivity of the LIGO network to the collection consisting of 
the five "most promising" millisecond pulsars, assuming 
they all have the same cllipticity. Our conclusions are 
summarized in Sec. [3 

II. NOTATION 

Let us consider an N-detector network, with output 
x a (t), a = l,...,N. (In cases where a single instrument 
outputs k independent data streams-e.g., a spherical bar 
detector that encodes for two GW polarizations in five 
data streams-we regard these formally as the outputs of k 
different detectors.) For simplicity, we begin by assuming 
that the detector noise is both stationary and Gaussian. 
We allow, however, for the possibility that the noises are 
correlated. Then we have 

< n a (f) n"(/'r >= \Kf - f')Sf(f) ■ (2.1) 

Here tildes denote Fourier transforms, according to the 
convention that 

/oo 
e 2 * ift x(t)dt, (2.2) 
-oo 

and "(■••)" denotes "expectation value". We note that 
S^f '(f) is the single-sided noise spectral density, which 
is also the convention followed in JKS. (If we were using 
the double-sided convention, the factor J on the rhs of in 
Eq. I|2.1[) would be replaced by 1.) 

The Gaussian random process n(f) determines a nat- 
ural inner product (. . . | . . .) on the space of functions 



x(t) 

( X \y) = mJ dfx<*(f)*[S^(f)] a/3 f(f), (2.3) 

where [S^ 1 (f)] aj3 Sf, 7 (/) = ^2 an d where 3? means "the 
real part of" . Here and below, to reduce index clutter, we 
sometimes represent a signal vector, having one compo- 
nent for each detector, by simply using boldface without 
an index; e.g., x(i) instead of x a {t). The inner product 
Eq. (|2.3|) is such that the probability distribution func- 
tion (pdf) for the noise n(i) takes the form 

pdf[n] = Me'^' 2 , (2.4) 

where N is a normalization constant. It follows that the 
expectation value of the product (x | n) (y | n) , over many 
realizations of the noise, is simply given by 

((x|n)(y|n)) = (x|y). (2.5) 



III. .^-STATISTIC FOR A DETECTOR 
NETWORK 

Given gravitational-wave data from a single detector, 
the .F- statistic developed by JKS is the optimal frequen- 
tist statistic for the detection of GWs from a single known 
NS in that single data stream. This section answers the 
question: if we have data from a network of detectors 
(possibly including bars as well as interferometers) how 
does one combine the different data streams to produce 
the optimal detection statistic for the entire network? 1 



A. jF-statistic for Single Source and Multiple 
Detectors, all with Time-Invariant Noise Curves 

Consider the search for nearly periodic GWs from a 
single source with known position and known (possibly 
time-varying) frequency, e.g., PSR 1937+21. The GW 
signal is characterized by four unknowns: an overall am- 
plitude Aq (equivalent to the combination /losin^sin 2 ^ in 
the notation of JKS), two angles i and ip that charac- 
terize the waves' polarization (equivalent to determining 
the direction of the NS's spin axis), and an overall phase 
$o- The GW signal h a (t) depends nonlinearly on l, ip, $ , 



JKS briefly consider this question and sketch a claimed answer, 
in §4 of their paper , but their answer is quite wrong. In par- 
ticular, they claim that the appropriate jF-statistic for a network 
with N detectors follows a \ 2 distribution with AN degrees of 
freedom, but we shall see below that the right number of degrees 
of freedom is just 4-the same as for the single-detector case. This 
is because there are still just four unknowns in the problem: the 
amplitude and phase of each of the two GW polarizations. 
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but, crucially, one can make a simple change of variablcs- 
to (A 1 , A 2 , A 3 , A 4 ) -such that dependence of h a (t) is linear 
in these new variables: 



sin 2 $(t), and cos$(t)sin$(t) by their time-averages: 



(3.1) 



where the four basis waveforms h" (t) are defined by 

h?(t) = F£(i)cos$ Q (t), hj(t) = F^(t)cos$ Q (t), 

fc°(t) = F*(t)sm$ a (t), hj(t) = F" (t)sin$ Q (t) . (3.2) 

Here $(t) is the waveform phase at the detector: 



$(t) « 2^ / / fltu (i')dt' 



(3.3) 



where f gw (t') is the measured GW frequency at the de- 
tector at time t'. The measured frequency includes the 
Doppler effect from the detector's motion relative to the 
source, as well as Einstein and Shapiro delays associ- 
ated with the Earth's orbit around the Sun. When the 
GW pulsar is in a binary, then f gw (t') also includes the 
Rocmcr, Einstein, and Shapiro delays associated with 
that binary orbit. We emphasize that the known-pulsar 
searches described here do not require the GW pulsar 
be isolated, but just that there exist an accurate timing 
model for the emitted waves. The F?(t) and F"(t) terms 
in Eq. 1)3. 2J) are the beam-pattern functions giving the re- 
sponse of the a th detector to the + and x polarizations, 
respectively. We note that the exact form of F"(t) and 
F"(i) depends on one's convention for decomposing the 
waveform into "plus" and "cross" polarizations; a one- 
parameter family of choices is possible, corresponding to 
the freedom to rotate the axes around the line of sight. 
JKS follow the conventions of Bonazzola & Gourgoul- 
hon 

A further word on our index notation: as above, we use 
Greek indices from the beginning of the alphabet (a, /?, 7) 
to indicate the various detectors in the network; we use 
Latin letters from the beginning of the alphabet (a, b, c) 
to indicate the four independent waveform components 
from a single NS (emitting at a single frequency); and we 
use Latin letters from the middle of the alphabet k) 
to label different NSs. As above, we sometimes remove 
the Greek index and instead represent the vector in bold- 
face: h (t) instead of h™(t). Finally, we use the capital 
Latin letter "J" to label different time intervals (always 
intervals over which the noise spectral density (/) can 
be safely approximated as constant); 

Next we define the 4x4 matrix T a b by 



_ , dh dh . , 



(3.4) 



Because both the observation time and 1 day (the 
timescale on which the F? x (t) vary) are vastly 
larger than the period of the sought-for GWs (typ- 



cos 2 $(t), sin 2 $(t) 
we have 



while cos$(t)sin$(t) — >■ 0. Then 



ri2 « E(^U-)) Q/3 /^(*)^(^ 

a, {3 

T 22 « Y.( S h 1 (f gW )) a J^(t)F^t)dt; (3.5) 

a, 

additionally, r 33 sa Ln, r 34 ps F i2 , r 44 w r 22 , and r i3 ps 
F14 pa F 23 sa F 24 sa 0. 

The best-fit values of A a satisfy 

^(x-^A b h b |.x-E ACh c) = (3.6) 



implying 



A a =^(F- 1 r fc (x|h b ) ; 



(3.7) 



and our optimal statistic 2T is then just twice the log of 
the likelihood ratio: 

2T = (x|x) - (x - E A"h b I x - J2 A c h c ) 

b c 

= E(r- 1 ) Qd (x|h a )(x|h d ). (3.8) 



Therefore using IT as one's detection statistic satisfies 
the Neyman-Pcarson criterion for an optimum test: it 
minimizes the false dismissal (FD) rate for any given false 
alarm (FA) rate. 

Writing x = n + h, and plugging into Eq. (|3.8() . we find 



(2T) = 4 + (h I h) , 



(3.9) 



where we have used Eq. 1)2. 5JI and the fact that 
< (h I n) > = 0. More generally, it is easy to show 
that y = 2!F follows a \ 2 distribution with 4 degrees of 
freedom (d.o.f) and non-centrality parameter p 2 = (h|h): 



P(y) = x 2 (yH;p 2 ) 



(3.10) 



ically 10" 



10- 



s), we can replace cos $(t), 



B. jF-statistic for a Single Source and Multiple 
Detectors with Time- Varying Noise Curves 

It is trivial to generalize the above results to a net- 
work of detectors with time- varying noise curves. Di- 
vide the total observation time into segments that are 
short enough that all noise correlation functions S^f 
can be approximated as constant during each segment. 
(We assume the segments are still much longer than the 
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GW period.) Let there be p such segments in all. Let 
the beginning and end points of these time intervals be 
(to,tx,--- ,t p ). (In this scheme, we can formally rep- 
resent gaps in the output of one or more detectors by 
intervals where some of the components S% go to in- 
finity.) While our signals come from N detectors with 
time-varying noise curves, we can formally regard them 
as coming from p N detectors, each with stationary noise 
(but such that only TV detectors are turned "on" at any 
instant; when TV more turn on, the previous TV turn 
off). But we know how to construct the JF-statistic for 
pN detectors with stationary noise characteristics, from 
the previous subsection. (Nothing in that subsection re- 
quired all the detectors to be "on" simultaneously.) The 
noise spectral density coefficients (/) are now labelled 
by time interval J: S? j(f)- Then Ln becomes 

T n ps EE(*))„J F?(t)FP(t)dt 

a,(3=U=l Jtj-i 

-> £ / P {S^(f gw ,t)) ap F^(t)F^t)dt (3.11) 

where we have made the notational shift S 1 ^ j(f gw (t)) — > 

S^(fg W (t),t); i.e., we have replaced the discrete label 
"J" by the continuous label "t". (In practice, the noise 
spectral density at any instant is estimated from the data 
itself, e.g., by use of a running mean.) 
Similarly, 

N f t P 

« E / {Sh\h^),t)) a0 n^) F0 At)dt 

N r l v 

r 22 w Y, I (S^(f gw {t),t)) ap F«{t)FP(t)dt (3.12) 

and again r 33 ps Tn, r 34 ps T 12 , r 44 ps r 22 , while T 13 ss 
Ti4 ps L 23 ps T 2 4 PS 0. 
If we define 

(x i h a ) = y r (Sh\f aw t)) af3 x a {t)hi (t)dt , 

a,/3 *° 

(3.13) 

then Eq. I|3.8|) remains the correct expression for 2JF, and 
Eq. H3.10fl remains its correct distribution function, with 
y = 2T and p 2 = (h|h). That is, given our notation, the 
expression for the multi-detector ^-statistic is the same 
as for the single-detector case. 



C. jF-stastic for a Single Source and N Detectors 
with Uncorrelated Noises 

The expressions simplify somewhat in the (common) 
case where the noises from different detectors are uncor- 
related: S^{f,t) = S%(f,t)5 a P. Then the inner product 



(x | y) is given by 



(x|y) 



' p x a (t)y a (t)dt 

'to 

We define A, B, C by 
As^xlhO, B=(h 2 |h 2 ), 



(3.14) 



C 



hi|h 2 ). (3.15) 



(Note that the A,B,C terms defined here are, in the 
single-detector case, larger than the A,B,C terms in JKS 
by a factor of the observation time To.) 

Then Ln = \A, L 22 = \B, and T 12 = \C. So L _1 
takes the form 



/ B —C 




(3.16) 



where D = AB — C . Thus we arrive at 



B{(x|h 1 )(x|h 1 ) + (x|h 3 )(x|h 3 )} 

+ A{(x|h 2 )(x|h 2 ) + (x|h 4 )(x|h 4 )} 
- 2C{(x|h 1 )(x|h 2 ) + (x|h 3 )(x|h 4 )}|3.17) 

As a check, consider the case of N identical, nearby de- 
tectors (assumed to have uncorrelated noises). Then A, 
B and C all scale like N, while D oc TV -2 . In the absence 
of a GW signal, the only terms in the (implied) double 
sum over a, (3 in (|3.17() that contribute, on average, are 
those with (3 = a. Thus terms like (x|hi)(x|hi) scale 
like N in the absence of a true GW, and so < 2T > re- 
mains invariant (always equalling 4) under changes of TV 
when there is no true signal. However when there is a 
true signal, then terms like (x|hi)(x|hi) scale like N 2 , so 
the non-centrality parameter p 2 of the distribution scales 
like TV-just as one would expect. 

Eq. H3.17f) can be re-written more compactly if we use 
complexified variables, as done in JKS. Defining 

2F a = (x|hx - ih 3 ) , 2F b = (x|h 2 - ih 4 ) , (3.18) 

Eq. (|3.17|) becomes 



IT = - [B\F a \ 2 + A\F b \ 2 - 2Cn(F a F b * 



(3.19) 



IV. ^-STATISTIC FOR MULTIPLE SOURCES 

In this section we consider a search for a collection of 
M nearly periodic GW sources, all with known positions 
and frequencies. In this case, the signal h a (t) depends 
linearly on 4M unknown parameters. Assuming that the 
M different GW frequencies /, (i = 1, • • • , M) are all 
sufficiently well separated that the detector noises are 
uncorrelated (i.e.,< n a (fi)h^{fj)* >= for i ^ j), then 
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a trivial repetition of the arguments in §111 shows that the 
optimum statistic (for either the single-detector or the 
multi-detector case) 2T is simply the sum of the optimal 
statistics for the individual sources: 

2T = J2 2J7 i- (4-1) 

i 

It also easy to show that y = 2T follows a % 2 distribution 
with AM degrees of freedom: 

P(y) = x 2 (yHM; P 2 )- (4.2) 

where the non-centrality parameter p 2 = p 2 . 

There are currently ~ 100 known millisecond (ms) pul- 
sars 2 (defined as pulsars with period P < 10 msec), of 
which ~ 60 are in binaries. We can of course consider 
any subset of these as some collection, sum their indi- 
vidual JF-statistics (derived from existing GW data) as 
in Eq. (|4.1|) . and test whether or not the collection has 
been detected. But when is such a strategy likely to 
be advantageous, and for which subsets? It seems that 
there are at least two interesting applications of this idea. 
First, one might hope that the nearest ~ 5 — 50 (say) ms 
pulsars, searched for as a collection, might be more de- 
tectable than any individual member. If this were the 
case, a multi-source search might hasten the first dis- 
covery of GWs from rotating neutron stars. We shall see 
below, however, that it is highly unlikely that a collection 
of more than a few (~ 2 — 5) of the brightest (in GWs) ms 
pulsars is more detectable than the very brightest source 
alone. While it is reasonably likely that the brightest few 
sources, taken together, are more detectable than the sin- 
gle brightest one-and we give a realistic example of this 
in Sec.IV.D-this will certainly not be the case for the 
brightest 20 or 50 sources. If there are too many sources, 
the strongest ones are effectively diluted by mixing them 
with the weaker ones, in the multi-source ^-"-statistic. 

To understand the second interesting application of 
multi-source searches, imagine a day when GWs have al- 
ready been detected from the few brightest, closest GW 
pulsars (all at distances of ~ 0.1 — 0.5 kpc, say), but 
when the GW pulsars in the range d > 0.5 kpc are 
still too faint to be detected. In that situation, it could 
make sense to take as a collection all (or some promising- 
looking fraction of) the ms pulsars in some annulus— say 
those in the range 0.5 < d < 1.0 kpc. This might al- 
low one to measure the average GW strength of those 
more distant sources, even if no single one of them could 
be positively detected in GWs, and therefore to begin 
to make interesting statistical statements based on this 
larger sample. 

We investigate the likely advantages of multi-source 
searches in the next five subsections. First, in Sec. IV. A, 



2 However ~ 60% of these are in globular clusters, at distances of 
several kpc, and so are roughly an order of magnitude further 
away than the closest known sources. 



we consider a collection of M sources, for large M, and 
ask: when does adding one more source to the collec- 
tion increase that collection's overall detectability? In 
Sec. IV. B we re-derive the distribution function of signal- 
to-noise-squared for any spatially uniform population of 
sources. These results from IV. A and IV. B are utilized 
in IV. C, where we show that for a uniform planar dis- 
tribution of GW pulsars (representing a somewhat ide- 
alized version of the population in our neighborhood of 
the Galactic disk), one might reasonably expect the few 
brightest sources, taken together, to be more detectable 
than the very brightest one. However this is hardly guar- 
anteed, and any advantages of a multi-source search are 
likely to be small in this case. This is illustrated in IV. D, 
where we consider a fairly realistic example based on the 
closet known ms pulsars. In IV.E we consider searching 
collectively for more numerous, more distant GW pulsars, 
after the nearest, brightest ones have been detected, and 
the advantages of multi-source searching are shown to be 
much greater in that case. 

A. The large-M case 

Here we compare the sensitivities of a single-source 
search and a search for a collection with M members, 
when M is much larger than one. For the single- 
source search, the threshold value of T that gives 
a 1% FA rate is given by 2T th (= yth) = 13.277 
(i.e.,J^ 277 x 2 (y|4)rfy = 0.01). To be detectable with 
FD rate > 50%, the signal strength must be at least 
p 2 = 10.234 (i.e.,/~ 277 x 2 (y|4; 10.234)% = 0.50). 

By comparison, when M is large, the x 2 distribution 
with AM d.o.f. is well-approximated by a Gaussian. Let 
y = 2T, and let p 2 tot = ^Li Pi- Thcn 

P(y) = X 2 (y\4M;p 2 ot ) « (8tt M)- 1 /2 e -(»-<»» 2 /(8A/) 

(4.3) 

where < y >= AM + pf ot . The thrcshhold value y t h such 

that f°° = 0.01 is then 

J yth 

y th ~ 4M + 4.652 VM (large M) . (4.4) 

(Note that the approximate threshold value that one ob- 
tains by inserting M = 1 into Eq. (|4.4(l is only 8.652, 
which is considerably less than the actual threshold yth = 
13.277 for the M — 1 case. Clearly, this is because the x 2 
distribution with only 4 d.o.f. has a substantial tail-i.e., 
is more skewed to the right than the higher- M distribu- 
tions.) 

When will a collection be more detectable than its sin- 
gle brightest member? To answer this, let us order the 
pulsars in the sample such that 

Pi > P2 > ■ ■ ■ > PM (4.5) 

Let Ti be the integration time necessary to detect the 
brightest source, and let T co u be the integration time 
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required to detect the M-member collection. For large 
M, the ratio of these two times is 



where 



TWTi = 0.455M 1 / 2 



A. 

n 2 

P tot 



0.455M~ 1/2 -^- , (4.6) 

r ave 



where p 2 ave = p\ ot /M . As an extreme example, if M = 
25 and all members have the same strength (pi — p2 ~ 
■ ■ • = P25), then Tcoii/Ti = 1/11.0. (More realistic cases 
will be considered in the next two subsections.) More 
generally, we say that a collection is more detectable than 
its brightest member if T co u/T\ < 1. 

How many pulsars should one include in the sample, 
assuming the goal is to hasten its detection? Imagine 
that pulsars 1, 2, • • • , M — 1 are included in the sample, 
and we want to decide whether to include pulsar M. By 
Eq. I|4.6[). the change AT co u in the time required to con- 
fidently detect the collection is 



AT, 



coll 



T, 



coll 



M-Ho.b-plJpl^ 



(4.7) 



Thus it is advantageous to increase the sample size (be- 
cause AT coU < 0) iff p 2 M > 0.5pl ve . 

Of course, a priori both p 2 M and p 2 ve are unknown; 
nevertheless one can use both some general statistical 
arguments and the measured parameters of nearby mil- 
lisecond pulsars to make a reasonably informed choice. 
We shall illustrate this in the next two subsections. 



B. Distribution of p for Galactic GW pulsars 

What is the distribution function of p 2 for the GW 
pulsars in the Galactic disk, within a few kpc of us? We 
can get quite far in answering this question, based on 
quite general considerations. 

Let r represent a pulsar's distance from the Earth. 
For simplicity, we shall consider two different spatial dis- 
tributions: a uniform (i.e., homogeneous and isotropic) 
three-dimensional distribution and a uniform planar dis- 
tribution. (These roughly represent the pulsar distribu- 
tions at distances r < 300 pc and 300 pc < r < 5kpc, 
respectively.) Let <j(r, /, A, oti) represent the probabil- 
ity density of GW pulsars in parameter space. Here / is 
again the NS's gravitational- wave frequency, A represents 
the signal's source's intrinsic amplitude (proportional to 



Eqw I '/ 2 , where Eqw is the source's GW luminosity), 
and the a, are the relevant angles in the problem. For a 
3-D distribution there are 4 such angles: two for the NS's 
angular location on the sky and two for the direction of 
its spin. For the planar (2-D) distribution there are only 
3 relevant angles, since one angle suffices for the sky loca- 
tion. For cither uniform distribution, the r— dependence 
can clearly be factored out: 



er(r, /, A, 014) = F(r)a(f, A, on) . 



(4.8) 



F(r) 



4irr 2 
2nHr 



for 3 
for 2 



D, 
D. 



(4.9) 



Here H sa 600 pc is the thickness of the Galactic disk. 

The source's signal-to- noise-squared, p 2 , can clearly be 
written in the following form: 



p 2 = A 2 r- 2 X(f, ai ) 



(4.10) 



where A(/, a,) is some function of / and the source's 
angular parameters. For a single detector with time- 
invariant noise characteristics, the /—dependence can 
also be factored out of A(/, aij: 



\(f, ai ) = X( ai )/S h (f); 



(4.11) 



however this factorization of the f-dependence is not nec- 
essary for our argument. 

For notational convenience, we again define y = p 2 . 
We now change variables: (r, /, A,ai) — > (y, /, A,ai). 
The density function a on the new variables is 



<7(y,f,A,ati) = <r{r, f,A,Oi) 



9(r, f,A,on) 



d(y,f,A,an) 



(4.12) 



where the second term on the rhs of Eq. Ij4.12|) is the 
Jacobian of the transformation. It is easy to check that 
this Jacobian factor is just 



d(r, f,A,ai) 



d(y,f,A,ai) 



1 



y-3/2^1/2 



(4.13) 



Combining Eqs. |@[HJ|, (f4"T!i|t and fiT^ . we there- 

fore have 

a(y,f,A, ai ) = &(f,A, ai )x { Jjla^f 7 * ' ^l^) 14 ) 

Integrating Eq. (|4.14|l over the variables (/, A, 014), we 
arrive at the density function for y alone: 



n 3 y 
n 2 y~ 



-5/2 



(3-D) 
(2-D), 



(4.15) 
(4.16) 



for some constants 713 and n%. 

We emphasize that no assumption about the distribu- 
tion of GW pulsars in f and A went into this result. All 
that was required was spatial uniformity-that the nearby 
pulsars are drawn from the same distribution as the more 
distant ones, and that the total number within some ra- 
dius scales as r to some power. Indeed, the same scaling 
applies to any source-type having a spatially uniform dis- 
tribution in Euclidean space; e.g., to the extent that one 
can ignore cosmological effects, the scaling in Eq. (I4.15f) 
also applies to detections of binary inspirals. (To appre- 
ciate why this is at least a bit remarkable, consider trying 
to estimate the density function cr(y) for all the pulsars in 
some globular cluster (say, 47 Tuc). In this case, all the 
pulsars arc effectively at the same distance, but we would 
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need to somehow estimate the distribution of GW pulsars 
in / and A, and then to fold in the detector's noise curve, 
in order to estimate a(y) for that cluster.) Of course, the 
scaling laws Eqs. (|4.15|l is well known in other areas of 
astronomy, where it is the basis for the ubiquitous log 
N-log S test of source strength distributions. 



C. Implications of o(y) for collection searches 

We can now return to the question of when a large 
sample of the brightest (in GWs) pulsars, taken together, 
might be more detectable than the single brightest mem- 
ber. From Eq. I|4.6[) . the "figure of merit" that character- 
izes the detectability of the collection is M~ x / 2 53 i=1 y 1 '. 
In the next two subsections we show how this quantity 
varies with collection size for spherical and planar distri- 
butions, respectively 



1. Uniform 3-D pulsar distribution 

The spherically symmetric case is the less interesting 
one, from a practical standpoint, since there only a few 
known ms pulsars with ~ 300 pc of the Earth. Neverthe- 
less we begin with this case since it is somewhat simpler 
and illustrates our general line of reasoning. 

Imagine that we have included in our collection all the 
brightest (in GWs) ms pulsars, down to some lower limit 
Vl . Then M pa /~ n 3 y~ 5 ' 2 and 



M 



Xy « / n 3 y-^dy, 



1/4 



(4.17) 



(4.18) 



This is a strictly increasing function of yi (albeit that it 
increases rather slowly). But increasing yi means shrink- 
ing the collection. Thus for a uniform 3-D distribution, it 
would be unlikely that a large collection of the brightest 
sources would be detectable before the single brightest 
member was detected. 



version of this is clearly 



M 
i=l 





-1/2 


/ 77,2 y~ 2 dy 




- Jyi 


Jyi 



I n 2 y~ 1 dy 

J 01, 



n\ /2 \n{y n / yi ) 
(n 22 /u) 1/2 



yr 1 



-1/2 



-hi(x) 



(4.19) 
(4.20) 



where in the last line we introduced the dimensionless 
ratio x = t/i/t/u < 1- We gain some insight into Eq. (|4.20|) 
if we re-express 772 in terms of y m ax, which we define to 
be the ?/- value of the strongest Galactic source. Let ymax 
represent the median value of y m ax, for our distribution 
function Eq. (|4.1ti|) . Then y max is given implicitly by 



77 2 y 2 dy = 0.5 . 



(4.21) 



(since then there is a 50% chance of finding a stronger 
source than y m ax), so y max = In-j. The actual value of 
ymax for our Galaxy is therefore y ma x = 2/3772, where 
we expect (3 is some number of order one. The rhs in 
Eq. H4.20f> can therefore be written as 



ymaxyu 



1/2 



-In (x) 



(4.22) 



1/2 

Next we consider the function /(x) = — lnx/(x _1 — l) , 
which is displayed in Fig. 1. 




2. Uniform 2-D pulsar distribution 

We turn now to the planar case. We begin by consid- 
ering the detectability of all GW pulsars with p 2 in the 
interval y± < p 2 < y u , (so yi and y u are the lower and up- 
per limits of the interval, respectively). Again, assuming 
the number of sources in the interval is large, the appro- 
priate figure of merit, characterizing the detectability of 
the whole collection, is M~ 1 / 2 'Y^L 1 y % ■ The continuous 



FIG. 1: Plot of function f(x) = -ln^)/^- 1 - l) 1/2 , which 
displays a maximum at x ~ 0.203. 

Note that it has a maximum at x ps 0.203, where 
/(x) pa 0.805. Thus an optimized source collection has 
a ratio of weakest-to-brightest source (in terms of their 
signal-to-noise-squarcd) of ~ 1/5. However the maxi- 
mum in /(x) is rather broad; at a brightness ratio of 
20 (x = 0.05), /(x) has decreased only ~ 15%, to 0.687. 
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Assuming the collection includes all the brightest sources 
down to some limiting brightest y 1 , the number of sources 
in the collection is 

M ~ n 2 / yi =y max /(2 yi ) (4.23) 

= 2.b[3- 1 ( Vmax,Vl ) . (4.24) 
5 

Thus including all sources down to strength y± while also 
optimizing ij\/y m ax leads to a rather small value of M, 
which is then somewhat outside the range of validity 
of the Gaussian approximation that led to our "figure 
of merit" il/" 1 / 2 YliLi V % 'i nevertheless it is clear that 
the "most-detectable" collection will have at most a few 
members. 

To estimate T co u/T\ (the time to detect this collec- 
tion divided by the time to detect the single brightest 
source), it would seem we need to evaluate the integral 
f n-22/ -1 (i-C, the continuous version of y l ), which 
is logarithmically divergent. Physically, though, it seems 
sensible to simply cut off the upper end of the integral at 
some y cu t of order y m ax- he., we cut off the integral at 
the y- value of the brightest source. 3 

Thus, setting y u equal to y m ax in Eq. (|4.22|) and plug- 
ging the result into Eq. Q4.6[l. we estimate 

» 0.455(2/3) 1 / 2 //(0.203) ss 0.80/3 1/2 . (4.25) 

Again, the use of Eq. I|4.6[) in deriving Eq. I|4.25[l is 
strictly valid only for large M; nevertheless the basic 
moral is clear: T co u/Ti is of order unity, and whether in 
actual experience it is greater or less than one depends 
strongly on (3, i.e., on whether the strongest source is 
stronger or weaker than one would expect, based on the 
source distribution function. 



D. Example: The Best Candidates among Known 
Millisecond Pulsars 

We next consider a potentially realistic example: a col- 
lection drawn from the population of known millisecond 
pulsars. We attempt to construct the "most detectable" 
collection from these. Of course, we do not know their 
actual GW strengths, so for this exercise we will esti- 
mate their strengths by assuming that they all have the 
same non-axisymmetry ie = I xx — I yy (where the NS 
is assumed to be spinning about its z-axis). This non- 
axisymmetry might be generated, e.g., by lateral varia- 
tions in the crustal composition or strong toroidal mag- 
netic fields in the NS interiors [l3[ . 



3 Of course, our planar approximation breaks down at r < H/2, 
and this "switchover" from an effective 2-Dtoa3-D distri- 
bution at short distances would obviate the need for an artificial 
cutoff in a more realistic treatment of this problem. 



For each of the millisecond pulsars, we estimate p 2 as 
follows. First, we estimate ho at the Earth from the 
pulsar's measured spin and the best available estimate of 
its distance r [l^, using 

h Q = 4ir 2 (G/c 4 )Ief 2 r- 1 (4.26) 

where here we will assume the GW frequency / is exactly 
twice the pulsar's measured spin frequency, v. Then we 
estimate p 2 using 

where To is some fiducial observation time and K{oii) 
is factor that depends on the sky-location and spin- 
oricntation of the source. The spin-orientations of the 
millisecond sources are poorly constrained, so for sim- 
plicity, in our estimates, we will simply replace K by its 
average value (over all angles). For Sh(f), we use the 
values for the Advanced LIGO noise curve, as generated 
by the Bench software package [l4|. (Eq. (|4.27|l is for a 
single detector; if one optimally combined the outputs of 
LIGO's LI, HI and H2 interferometers, then p 2 should 
be approximately 2.25 times greater than for either LI 
or HI alone.) 

Given the above inputs, we find that there are 5 ms 
pulsars that stand out as the best candidates for detec- 
tion by the Advanced LIGO Interferometers. They are 
PSRs J0437-4715, J0030+0451, J2124-3358, J1744-1134, 
and J1024-0719. These 5 are also the closest known 
ms pulsars. And pulsar 1 (PSR J0437-4715), which at 
d = 0.14 kpc is the closest of all the known ms pul- 
sars, is estimated to be the strongest GW source. Rel- 
ative to pulsar 1, the GW strengths of the other four 
sources are given by: (/O2//O1) 2 = 0.38, {pz/pi) 2 = 0.32, 
(p 4 /pi) 2 - 0.17, and (p 5 / Pl ) 2 = 0.16. 

Assuming the above estimates of p 2 for the five 
best candidates were correct, what would be T co u/T\ 
(the ratio of the integration times necessary to de- 
tect the 5-member collection and the brightest individ- 
ual source)? For our 5-member sample, the threshold 
value for detection with 1% FA rate is yth = 37.57 
(i.e.J™ 57 x 2 (y\20)dy = 0.01). To be detectable with 
FD rate > 50%, the signal strength must be at least 
pl t = 18.45 (i.e.,/ 3 ~ 57 x 2 (y|4; 18A5)dy = 0.50). Thus 

| = W)(|^)=0.87. (4.28) 

E.g., if it took two years to confidently detect the 
strongest source, the 3-member ensemble would be de- 
tectable in about 21 months. (Note that in deriving 
Eq. (|4.28(l we have used the actual x 2 distribution with 
20 d.o.f., not the Gaussian approximation to it.) 

Should we add a sixth pulsar to the sample? From 
the analysis in the previous subsection, this would be 
advantageous if p§/(£i=i Pi) > 0.10. But we estimate 
that the sixth most detectable pulsar is J1012+5307, with 
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(pe/pi) 2 = 0.07, so we restrict the sample to the most 
promising five. [Indeed, if we had restricted the sample to 
only the most promising 3 pulsars, we would coinciden- 
tally have arrived at the same estimate for T co ujT\. For 
the 3-member case, (Ptot/Pi) = 1-70, while the thresh- 
old for detection with 1% FA rate is yth = 26.22, and 
the signal strength must be p 2 ot > 15.13 to be de- 
tectable with a FD rate > 50%. Thus we would estimate 
Tcoii/Tx = 15.13/(2.03 * 10.23) = 0.87, the same as for 
the 5-member collection. However we highlighted the 5- 
member result since that one is clearly somewhat more 
robust against deviations of the actual source strengths 
away from our fiducial estimates.] 

Clearly, since the actual orientations of the ms pulsars 
are unknown and the distances are known only to within 
a factor ~ 2, the above estimate merely gives a rough 
indication of the time-savings that a multi-source search 
might reasonably lead to. 

We also note that if we were to estimate source 
strengths by assuming that all ms pulsars are spin- 
ning down primarily due to GW emission, thus using 

ho = (^7) 1/2 ^ _1 instead of Eq. I|4.26|l . and then repeat 
the above analysis from that starting point, we would find 
that no subset of the known ms pulsars is more detectable 
than the single brightest source, PSR J0437-4715. This 
just highlights the fact that the ratio T co /;/Ti-and es- 
pecially whether that ratio is greater or less than one- 
depends rather sensitively on the relative strengths of the 
few brightest sources, which of course we will not know 
in advance. 



E. Search for a collection of weaker GW pulsars 

The last two subsections showed that a search for a col- 
lection of the very brightest GW pulsars may offer some 
advantages, compared to a search for the very brightest 
one, but any such advantages arc likely to be quite mod- 
est. We now turn to a case where the advantages of a 
whole-collection search are much more impressive. 

Consider some time in the future, when the few bright- 
est GW pulsars are presumed to have already been de- 
tected. These are presumably among the closest GW 
pulsars, while the GWs from their more distant cousins 
are still too weak (at the Earth) to be detected. Now once 
again consider collecting together all known ms pulsars 
in the range yi < y < y n . (Of course, again, one does not 
know precisely which these arc, but whatever lessons are 
learned from the brightest GW pulsars, combined with 
the known distances and spin rates of the remaining mil- 
lisecond pulsars, will probably allow one to make fairly 
educated guesses.) If all ms pulsars had the same intrinsic 
GW strength, the same frequency, and the same angular 
factor K{cti) (from Eq. I4.27|) . then clearly these pulsars 
would occupy a circular annulus in the disk. In fact, of 
course, these GW pulsars will not have the same intrinsic 
amplitude, frequency, or angular factor, but we still find 
it helpful, conceptually, to imagine the GW pulsars with 



?/i < y < y u as filling a roughly annular region. 

For an optimally chosen annular region (one that min- 
imizes the integration time required for positive detec- 
tion), what is the optimal value of x = yi/y n - We 
worked this out in Sec. IV. C. 2; the optimum selection has 
x ~ 1/5. How many sources are in this range? For a 
planar distribution, the answer is clearly 

M« r /5 n 2 y- 2 dy = 2(^), (4.29) 

Jy u Vu 

where we have used n 2 = y m ax/2- Similarly the ratio 
T C oii/T u (where T u is the integration time required to 
detect a GW pulsar whose p 2 equals y u ) is given by 

T coU /T n » 0.455i; u /[K2/u) 1/2 /^«-203)] (4.30) 
« 0.8 (J^) 1/2 . (4.31) 

yraax 

For example, if y u /y m ax = 1/10, then the optimal num- 
ber of sources for the "annulus" is M w 20, and the in- 
tegration time required to detect that whole collection 
of 20 GW pulsars is only 0.8/\/l0 « 0.25 as long as the 
time required to detect the brightest single member in 
that group. Again, such a detection would provide an 
estimate of the average p 2 for ms pulsars in that collec- 
tion, even though none could be detected individually in 
GWs. 



V. CONCLUSIONS 

The most sensitive GW detectors (currently the LIGO 
LI, HI, and H2 interferometers) have very similar sen- 
sitivities, and this is likely to remain the case for some 
years. In such a case, one can significantly increase the 
effective signal-to-noise of any source by optimally com- 
bining the data streams. Here we have derived the ap- 
propriate formulae for doing so, for GW pulsar searches. 
For N GW detectors with the same sensitivity, the ob- 
servation time Tdet required to detect any particular GW 
pulsar scales like ./V -1 , and so combining the data streams 
this way is clearly a useful strategy. 

However we remind the reader that our analysis of the 
multi-detector statistics has assumed the noise is Gaus- 
sian. In the more realistic case, one would one want to 
veto candidate detections that had a large .F-statistic but 
that did not sufficiently resemble actual GW pulsar sig- 
nals, e.g., because the relative sizes of the signal in the 
various detectors did not conform with expectations for 
any choice of parameters (A 1 , A 2 , A 3 , A 4 ). In particular, 
we imagine that a realtistic implementation would in- 
corporate some multi-detector version of the chi-squarc 
veto developed in Itoh et al. 01 1 however we have not 
considered this in any detail. 

We have also pointed out that one can search for collec- 
tions of pulsars, and that the optimal frequentist search 
for such collections simply adds up the ^"-statistics of 
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the individual members. We considered two cases in 
detail. We first asked whether the few brightest GW 
pulsars might be discovered, collectively, before the very 
brightest one. The answer turns out to depend rather 
sensitively on the relative strengths of the few bright- 
est sources, and so we can only equivocate: maybe yes, 
maybe no. But even if some collection turns out to be 
more detectable than the single brightest source, it is un- 
likely to "win" by much. However, after the few brightest 
GW pulsars have been discovered, searching for more dis- 
tant pulsars by summing their ^"-statistics should prove 
to be an effective strategy, allowing one to measure the 



average strength of many sources that are not individu- 
ally detectable. 
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